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1. Introduction 

Many models for physics beyond the Standard Model, notably Super symmetry, Little Higgs 
models, various models for extra dimensions and several versions of technicolor models, 
contain new strongly interacting particles with masses below or near the TeV scale. These 
particles, if they exist, will be copiously produced at the LHC, and could provide the first 
and most important sign of new physics. By analyzing their production and decay modes, 
we might be able in only a few years to determine many of the new physics properties, 
including large parts of the particle spectrum, spin structure and the possible existence of 
new stable particles that can act as dark matter candidates. 

The final state signatures characteristic of events where heavy colored particles are pro- 
duced typically include a large energy measured in the detector, missing transverse energy 
and several hard jets coming from their decays. Such very spectacular events can be effi- 
ciently simulated by Monte Carlo generators that implement new physics models, including 
general-purpose generators like Pythia Herwig Q and Sherpa ||, and matrix ele- 
ment generators such as CalcHEP/CompHEP j|, ||, MadGraph/MadEvent ||] and 
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Whizard Q. However, an additional difficulty in the simulation of particle production at 
hadron colliders is due to the presence of abundant QCD radiation, in particular initial- 
state radiation. This radiation, which is enhanced in the production of heavy and strongly 
interacting particles, can have important effects. It affects the event kinematics by giving a 
transverse boost to the heavy particle system, and can produce additional jets besides the 
jets originating from the decay of the heavy particles, thus complicating the reconstruction 
and indentification of the event. 

This additional jet production has traditionally been simulated using Parton Shower 
(PS) Monte Carlo programs such as Pythia and Herwig, which describe parton radiation 
as successive parton emissions using Markov chain techniques based on Sudakov form 
factors. This description is formally correct only in the limit of soft and collinear emissions, 
but has been shown to give a good description of much data also relatively far away from 
this limit. However, for the production of hard and widely separated QCD radiation jets, 
this description breaks down due to the lack of subleading terms and interference. For that 
case, it is necessary to use the full tree-level amplitudes for the heavy particle production 
plus additional hard partons. 

The Matrix Element (ME) description diverges as partons become soft or collinear, 
while the parton shower description breaks down when partons become hard and widely 
separated. In order to describe both these areas in phase space, the two approaches must 
be combined, without double counting or gaps between different parton multiplicities. An 
additional physical requirement is that such a procedure should give smooth distributions, 
and interpolate between the parton shower description in the soft and collinear limits and 
the matrix element description in the limit of hard and widely separated partons. Several 
algoritms have been proposed to achieve this, including the CKKW Lonnblad ||To|| 

and Mangano [11, 12] schemes. These different procedures are in substantial agreement 
and give consistent results at hadron colliders For Standard Model processes, in 

particular the production of jets and weak vector bosons, the matching schemes have been 
extensively used and compared to the available data from the Tevatron ||15[| . In addition, 
inclusive top quark pair production with matching has been compared to a next-to-leading 
order plus parton shower approach (MC@NLO) Jl6| in Ref. fl~2f| . 

That jet matching is necessary for Standard Model backgrounds is obvious, since the 
only way many Standard Model processes (such as weak vector boson production) can 
emulate the production of new heavy states charged under QCD is by taking into account 
high-pj_ radiated jets. It is less clear, however, that jet matching should be important in 
new physics production, for several reasons: The new particles are expected to decay to 
high-pj_ jets making them fairly insensitive to the effects of QCD radiation; furthermore 
the parton shower formalism should work better the higher the mass is for the produced 
particles, since then the region in which emissions can be considered "collinear" is increased. 
As we will demonstrate in this paper, this is only partly true, and there are many scenarios 
where jet matching turns out to be surprisingly important also in the production of heavy 
new physics QCD states. 

We will here use the production of squarks and gluinos within the framework of the 
Minimal Supersymmetric Standard Model, MSSM. The properties of initial state QCD 
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radiation are however quite insensitive to the precise type of particle produced, and mainly 
depend on the mass scale and production mechanism of the new particles. All our results 
here are therefore readily applicable to any type of new physics producing colored heavy 
states decaying to jets. 

The layout of the paper is as follows: In Sec. ^ we discuss the details of the matching 
schemes used, as well as how to solve the double counting problem that arises due to the 
decay of on-shell states in multiparton final states. In Sec. ^ we study the differences in 
QCD radiation between matched and unmatched simulations. The matching is found to 
reduce the uncertainties of the parton shower and hence to increase the predictivity of 
the simulations. This increased predictivity allows us, in Sec. |3], to study the effects of 
QCD radiation in the presence of jets from the decay of the new states. We analyse in 
detail a few examples where jet matching turns out to be crucial to get a qualitatively 
correct description of the new physics signatures. In Sec. || we take a look a the impact 
of using matched simulations in a more realistic experimental situation, where we include 
production of all relevant supersymmetric particles as well as Standard Model backgrounds. 
We draw our conclusion in Sec. ||. 



2. Technical challenges 

In this section we present the technical issues and the problems which can arise while merg- 
ing matrix elements involving heavy colored states with a parton shower. Our approaches 
and their implementations are completely general and can therefore be applied equally well 
to the SM as to any Beyond the SM (BSM) construction. To be concrete we will focus on 
top pair production and on pair production of strongly interacting SUSY particles, such 
as squarks and gluinos. In this study we have employed the MadGraph/MadEvent 
matrix element generator jl7], || interfaced to Pythia 4 [|]J for parton showering and 
hadronization. The MSSM implementation used was presented in Refs. |2(J [|. 

2.1 Matching schemes 

The goal of jet matching is to merge samples with different parton multiplicity obtained 
via matrix elements, correctly accounting for showering effects and avoiding double count- 
ing. The matching algorithms used in this study can be viewed as hybrids between the 



approaches currently employed by Sherpa [|| and Alpgen [21]. The phase space sepa- 
ration between the different multijet processes is achieved using the /jj_-measure p2]] . No 
analytic Sudakov reweighting of the events is performed, but instead showered events are 
rejected if they are not matched to the matrix element-level partons. This method allows 
the use of the well-tuned showering and hadronization implementations of Pythia while 
retaining the advantages of matrix element production. The matching implementations in 
MadGraph/MadEvent can be used both for Standard Model and new physics processes. 
The first matching scheme used in this study, the fcj_-jet MLM scheme, is the one 



used for MadGraph/MadEvent in Ref. Q. The final-state partons in an event are 



clustered according to the k±-jet algorithm to find the "equivalent parton shower history" 
of the event. Here, the Feynman diagram information from MadGraph is used to allow 
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only clusterings that correspond to diagrams existing in the generated matrix element. The 
smallest k± value is restricted to be above some cutoff scale ■ In order to closely mimic 
the behaviour of the parton shower, the k± value for each clustering vertex corresponding 
to a QCD emission is used as renormalization scale for a s in that vertex. As factorization 
scale, as well as renormalization scale for the central hard 2 —* 1 or 2 —* 2 process, the 
transverse mass = p\ + m 2 of the particle(s) produced in the central process is used. 
This event is then passed to Pythia for parton showering. After showering, but before 
hadronization and decays, the final-state partons are clustered into jets, again using the 
k± jet algorithm, with a cutoff scale Qmatch > Qcut • These jets are then compared to the 
original partons from the matrix element event. A jet is considered to be matched to the 
closest parton if the jet measure k± (parton, jet) is smaller than the cutoff Qmatch- The event 
is rejected unless each jet is matched to a parton, except for the highest multiplicity sample, 
where extra jets are allowed below the k± scale of the softest matrix element parton in the 
event, <2^ft est - This matching scheme can be used with both the old (vituality-ordered) 
and the new (p^-ordered) shower implementations of Pythia. 

In order to further study the systematics of the jet matching, a new matching scheme 
has been implemented, which we call the "shower k± scheme". In this scheme, which is 
here used for the first time, events are generated by MadGraph/MadEvent as described 
above, including the reweighting of a s . The event is then passed to Pythia and showered 
using the pj_-ordered showers. For the pj_-ordered showers, Pythia reports the scale of the 
first (hardest) emission in the shower, ^hardest- For events from lower-multiplicity samples, 
the event is rejected if Qhardest ^ s above the matching scale Q ma tch) while events from 
the highest multiplicity sample are rejected if Qhf rdest > Q^ft es t> the scale of the softest 
matrix element parton in the event. This matching scheme is simpler and yet effectively 
mimics the workings of the fej_-jet MLM scheme. However, it allows for the matching scale 
Qmatch to be set equal to the matrix element cutoff scale Q^f, and it more directly samples 
the Sudakov form factor used in the shower. Furthermore, the treatment of the highest 
multiplicity sample more closely mimics that used in the CKKW matching scheme. 

Since we here study heavy QCD particle production at hadron colliders, we expect 
final-state radiation to be highly suppressed as compared to initial-state radiation. The 
treatment of final-state radiation from the (undecayed) produced particles is therefore of 
somewhat secondary importance, and does not strongly affect any results of this study. 
We have here chosen not to differentiate between shower particles originating from final- 
state radiation and initial-state radiation in the jet clustering in the k±-jet MLM scheme, 
and correspondingly for the "shower k±" scheme, veto also events where the final-state 
radiation reports an emission scale above the matching scale. 

2.2 Double counting from resonant diagrams 

In the simulation of events where colored particles are produced that can decay into one 
another by emitting partons, a further problem arises, also related to double counting. 
Consider, for example, the contributions to gq + 1 jet coming from various subprocesses. 
Among these, gg — > gqq displays a peculiar behaviour. In Fig. [j] we show three representa- 
tive diagrams out of a total of sixteen. Diagram (a) is a "genuine" correction to the 2 — > 2 
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Figure 1: Representative Feynman diagrams contributing to gq+l jet via the subprocess gg — > gqq. 
Diagram (a) is a "genuine" QCD correction to the Born level diagram and correctly handled by the 
matching procedure. Diagrams (b) and (c) belong to the same gauge invariant class as (a), however 
they contain a possibly resonant gluino (squark) contribution, relevant when rrig > rriq {rrig < rriq). 
These processes are already taken into account in the gg and qq* channels, respectivetly, and have 
to be subtracted. 

Born amplitude and it is correctly handled by the matching procedure. Diagrams (b) and 
(c), however, contain possibly resonant gluino and squark propagators. When integrated 
over the phase space these diagrams give rise either to gg — ► gg with g — ► qq or to gg — ► qq* 
with q* — > gq, depending on the mass hierarchy. These contributions, however, are already 
taken into account by the Born level processes gg — ► gg or gg — > qq* plus the corresponding 
decays and have therefore to be properly subtracted to avoid double counting. In fact, this 
kind of issue is not specific to SUSY, but appears every time a given final state can be 
reached through different decay cascades. 1 Several solutions have been proposed in the 
literature, with various degrees of approximations, only a few of which can be employed in 
a Monte Carlo framework. For example, subtractions at zero width, typically implemented 



in NLO calculations, are not suitable to a Monte Carlo approach. In Ref. [25] these di- 
agrams were removed from the matrix element by hand. Even though this procedure in 
principle violates gauge invariance, it works very well in practice when the width of the 
resonances is sufficiently small, T/m <C 1. In order to deal with this kind of problem in 
full generality we have implemented and checked two independent types of solutions. 

In the first, resonant diagrams are removed as in Ref. p5| , but in an automatic way 
directly by MadGraph. This approach has the virtue of being very simple, though it is not 
gauge invariant and neglects the interference between the diagrams. The impact of these 
approximations therefore has to be carefully checked case by case. 

The second solution is an hybrid between a proper gauge-invariant subtraction and 
the diagram removal. It is based on a general algorithm implemented in MadEvent to 
record the information of the presence of intermediate resonant propagators on an event- 
by-event basis. Such information is needed by the showering program to ensure that the 
QCD evolution does not shift or even cancel the Breit-Wigner peaks. When only resonant 



A well-known example in the Standard Model is in the calculation of the strong corrections to tW 
production. In this case the process gg — * tWb contains a diagram with a resonant top which overlaps with 
tt production E| E3] 
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or only non-resonant diagrams are present, there is, in principle, no ambiguity and the 
information could be directly passed to the events. When both resonant and non-resonant 
diagrams are present such a unique assignment is not possible, due to interference effects. 
However, the interference terms are typically suppressed by T/m, and can therefore be 
neglected in the narrow width limit. In this case one can associate the presence of a 
propagator to an event on a statistical basis, using the relative size of the corresponding 
squared amplitudes. The above procedure exactly parallels that commonly employed to 
associate a color flow to an event [26]: The calculation of the matrix element is exact at all 
orders in the number of colors (as the amplitude-squared includes all diagrams and their 
interference), however the color flow assignement is approximated by ignoring interference 
terms which are proportional to 1/N%. Such propagator/diagram mapping is possible 
in MadGraph/MadEvent due to the fact that automatic integration over space space 
is already based on the amplitudes corresponding to single diagrams. The phase space 
integration is split into independent channels, each corresponding to a single Feynman 
diagram Events generated in a resonant channel are provided with the resonant 

particle information and also one of the color flows consistent with the Feynman diagram 
itself (in so doing, we avoid possible mismatches between the propagator structure and the 
color flow, which would lead to an inconsistent shower evolution). Having at our disposal 
the above procedure, it easy to get rid of the resonant contributions, by simply dropping the 
events with resonant propagators before passing them to the Pythia interface. The results 
presented and discussed in the following are obtained using the latter method. However, 
we have checked that the two methods are consistent. 



3. QCD radiation in ME/PS merged samples 

QCD radiation in production of heavy states at hadron colliders is expected to be mainly 
due to initial state radiation. 2 The main reason for this is that the phase space for initial 
state radiation is typically large and that collinear radiation from the final state particles 
is suppressed by their large mass. The main impact of QCD radiation on the production 
of heavy states is therefore twofold: A transverse boost of the heavy particle pair, and 
production of additional possibly hard jets radiated from the incoming parton lines. The 
details of the initial state radiation are expected to depend mostly on the mass of the 
produced particles, which sets the scale for the radiation, and on the type of initial state 
partons (whether the production is from valence quarks, quark-antiquark or gluon fusion). 

In this section we compare the QCD radiation pattern as obtained from a parton shower 
alone (Pythia) and from a matrix element plus parton shower (MadGraph/MadEvent 
+ Pythia) approach, for some key processes with heavy colored final states. To sim- 
plify the discussion and reduce the model dependence, we study specific heavy final states 
which are left undecayed, and give all distributions after parton showering but without 
hadronization. We compare the systematic uncertainties involved in parton showering 

2 While this is strictly speaking a non gauge-invariant statement, given that radiation is associated to 
color flows and not to particles, it has become common jargon to talk about initial and final state radiation. 
Such a separation is in fact meaningful only for the collinearly enhanced radiation. 
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by using the two different showering implementations of Pythia, the "old" (virtuality- 
ordered) and "new" (pj_-ordered) showers, with a range of shower parameters similar to 
that in |p5| . Similarly, for the matching we use the same shower implementations and pa- 
rameter variations, and furthermore use two different matching schemes, the fcj_-jet MLM 
scheme for virtuality-ordered showers and the "shower k± scheme" for p_|_-ordered showers, 



as discussed in Sec. 2.1 



For most of the studies below, the shower parameter that is varied is the starting 
scale of the shower. This is the most important parameter to determine the hardness of 
radiation allowed in parton shower emissions, and the default value has varied over the 
years. Originally, the default scale was set to the factorization scale /A of the process 
(typically m\ = p\ + rn 2 of the produced particles), inspired by the notion of (initial 
state) parton showers being the deconvolution of the DGLAP equation for the parton 
density functions. However, with the so-called Tune A p7| , based on W and Z boson 
production at the Tevatron, this was replaced by 4//|*. Later (from v. 6.319), this upper 
bound was increased, for certain types of processes, to s pp , the center of mass energy of 
the proton- (anti) proton collisions, in order to further improve the description of extra jet 
radiation in hadron collisions. Following the nomenclature of Ref. pq] , we will in the 
following call showers with a starting scale corresponding to \j? f "wimpy showers", and 
showers with the starting scale s pp "power showers". 

To study the impact of QCD scale choices, we also vary the factorization and renor- 
malization scales by a factor 1/2 and 2, in the matrix element as well as the parton showers. 

We consistently simulate X+ jets for X = gg,q^q^*\gq^ and tt, with up to two 
extra partons from the matrix element simulations. 

3.1 Variation of matching parameters 

Before starting the comparisons between unmatched (showered-only) and matched distri- 
butions, a series of "sanity checks" on the matched distributions were performed. These 
include the requirement of a smooth transition between the region of phase space described 
by the shower (below the matching scale Q ma tch) and the region described by matrix ele- 
ments (above Q ma tch); and stability of distributions as the matching scale is varied. Since 
both the matching implementations employed here rely on the Durham k± measure to 
achieve the separation of the phase space, the most revealing distributions to study their 
features are the differential jet rates defined according to the same measure. In particular, 
in the fc_i_-jet MLM scheme there is at parton level a sharp division in the jet rates between 
the shower and matrix element regions, making it very easy to see to which extent the 
transitions are smooth. For the shower k± scheme, as well as the cone jet MLM scheme 
implemented in Alpgen and the CKKW scheme implemented in Sherpa, the separation 
is less sharp, but the differential jet rates still tend to be the best variables to study the 
transition between parton showers and matrix elements. 

Our guidelines for the choice of the scales Q^f and Q ma tch are based on smoothness of 
distributions across the matching transition as well as efficiency. The higher the scale can 
be chosen, the higher will the proportion of lower-multiplicity events be in the combined 
sample. Since the computational effort is heavier for higher multiplicities, it is desirable to 
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choose the matching scale as high as possible. Therefore for each choice of particle type, 
particle mass, and shower, the matching scale is chosen to be close to the highest scale that 
still gives a stable matched cross section and smooth differential distributions, in particular 
differential jet rate distributions. As it turns out, the "new", pj_-ordered Pythia showers 
allow significantly higher choices of the matching scale than the "old", virtuality-ordered 
showers. The reason for this is that they give significantly harder emissions than the old 
showers, and therefore give distributions more similar to the matrix element distributions. 
In the following, our default scale choices are: 

• For pair production of SUSY particles with heavy mass (600 GeV and above): For 
virtuality-ordered showers, = 40 GeV and Q ma tch = 60 GeV, for pj_-ordered 
showers Q^E and Q match = iqo GeV. 

• For Standard Model ti production or light gluino par production: For virtuality- 
ordered showers, Q^f = 20 GeV and (Jmatch = 30 GeV, for p^-ordered showers Q^f 
and Qmatch = 100 GeV. 

In Fig. I we show the effects of varying the matching scale by a factor 1/2 and 2 
from the default values, for 607 GeV gluino pair production and Standard Model top pair 
production, for the new and old Pythia showers. We have here used the default starting 
scales for the respective Pythia showers, corresponding to the factorization scale \if for 
the pj_-ordered showers and for the virtuality-ordered showers. 

In the right-hand column, the curve for unmatched Pythia showers with default pa- 
rameters are shown, together with the pure matrix element prediction without any parton 
showering or matching and the matched curve. We see, as expected, that the matched 
curve smoothly interpolates between the unmatched Pythia curve below the matching 
scale, and the matrix element prediction for large scales. 

We also see how the variation in the matched curves is small as the matching scale is 
varied, and the curves above and below the matching scale variation limits are quite stable. 

The difference between the two matching schemes is visible in the plots, in the different 
behaviour of the parton multiplicity sample contributions. The lefthand column for each 
particle type shows the the fcj_-jet MLM matching scheme, with the contributions from 
the different parton multiplicity samples in grey. The matching scale cutoff is, in this 
scheme, done in the same variable that is plotted, the differential jet rate, and there is 
therefore a sharp cutoff between the 0- and 1-parton samples in DJR(0 — * 1), and between 
the 1- and 2-parton samples in DJR(1 — > 2), so that below the cutoff only the lower- 
multiplicity samples contribute and above the cutoff only the higher-multiplicity samples. 
In the middle column, the "shower k±" scheme is used (with the p^-ordered Pythia 
showers). This scheme cuts on the first emission of the parton shower rather than on the 
combined radiation of the whole shower, giving some smearing across the matching scale. 
This scheme therefore allows to use the same cut at matrix element level and matching 
level. The distributions for the p^-ordered showers have been double-checked using the 
k±-jet MLM matching method, with excellent agreement. 
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Figure 2: Differential jet rates for — > 1 and 1 — > 2 jets from QCD radiation for gg and ti 
production at the LHC. The first two columns show the distributions for the two types of Pythia 
showers, virtuality-ordered and pj_ -ordered, for three different choices for the matching scale. For 
the default choice, also the contributions from the separate multiplicity samples are shown. The 
colored arrow show the value of the three Q ma tch used for each kind of shower. The third column 
shows how the matched curve interpolates between the pure parton shower curve and the pure 
matrix element curve without parton showering. 



Interesting to notice is the differences in curve shapes depending on the choice of shower 
type. Below the matching scale, the shape of the curve is given completely by the shower, 
in particular for the — > 1 jet rate. Above the matching scale, however, the shape is 
mainly given by the matrix element. It is easy to see the reason for the different choices of 
matching scales for the different showers - the pj_-ordered shower gives significantly harder 
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Figure 3: p± spectrum for the first and second hardest radiated jet in (a) gg and (b) ti events. The 
right column in each group of plots shows the spread of Pythia predictions with different choices 
of the shower evolution variable (virtuality- and p^-ordered) and starting scale for the evolution 
(labeled as "wimpy" and "power" showers respectively). The left column presents the results 
obtained after matching in the same four shower scenarios. The grey band shows the systematic 
uncertainty associated with a variation of a factor of two of the renormalization and factorization 
scales with respect to their central values. The different curves have a normalization corresponding 
to their cross section. The gluino mass is here 607 GeV, while the top mass is set to 174 GeV. 



distributions than the virtuality-ordered shower, and is more similar to the matrix element 
curve, hence allowing a higher matching scale. 

3.2 Parameter dependence in matched and unmatched generation 

One of the advantages of the parton shower formalism, and also one of the arguments 
for using it, is that parton shower Monte Carlo generators have several parameters on 
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which the behaviour of the shower depends, that can be tuned to the data. While this is 
certainly an advantage in general, it also means that the parton shower lacks predictability 
at least for some observables or areas of the phase space. It is not always clear that a tune 
done for one type of initial state will be applicable to other initial states, or that a tune 
done for a particular mass of a pair-produced particle will be applicable for other masses. 

We first compare the matched and unmatched spectra of the first and second jet p± 
in Standard Model tt production events and SUSY gg production with nig = 607 GeV, for 
different shower parameters, Fig. |3[ Here, we look at only jets coming from QCD radiation, 
leaving the heavy particles undecayed, in order to specifically study the differences in treat- 
ment of radiation only. The different curves correspond to the different shower parameters 
settings: the "old", virtuality-ordered shower and the "new", pj_-ordered implementation, 
each with two different choices for the starting scale of the shower, [ip ("wimpy shower") 
and s pp ("power shower"). In order to study the intrinsic QCD uncertainty on the predic- 
tions, we also vary the factorization and renormalization scales for one of the parameter 
settings by a factor 2 up and down. Note that we here vary all scales together, both for 
the central process and for the parton shower /QCD radiation in the matching. 

Several interesting features can be noted from Fig. ^ 

• The spread in predictions for the parton shower is very large and strongly affects 
the shapes of the distributions. This uncertainty due to shower parameters is almost 
completely removed when matching is applied. 

• The region where the shower predictions start to diverge, and the rate of this diver- 
gence, is strongly correlated with the mass of the produced particles. This correlation 
is due to the choice of starting scale for the "wimpy" showers as the factorization 
scale, which is close to the mass of the produced particle. 

• The "power" shower curves consistently overshoot the matched curves, and hence 
give too hard predictions, while the "wimpy" showers give too soft distributions. 

• The uncertainty due to scale variations is considerable, but mainly affects the nor- 
malization and only to a small degree the shape of the curves. 

The present Pythia default for the virtuality-ordered showers is close to the curve for 
the "wimpy" pj_-ordered showers, which is also the default for the pj_-ordered showers. 

3.3 Impact of different inital states 

The impact of different initial states is shown in Fig. ||, where we compare radiation for pair 
production of gluinos, up- type squarks and up- type squark-antisquark. For ease of com- 
parison, the supersymmetric particle masses have all been set to 607 GeV. These different 
SUSY particles are produced by different initial states: gluino production is mainly through 
gluon fusion, squark-squark production from valence squarks exchanging a t-channel gluino, 
and squark-antisquark production from a combination of gluon fusion and quark-antiquark 
with t-channel gluio exchange. Besides the pj_ of the two hardest radiated jets, we also 
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Figure 4: an d P± distributions for radiated jets in gg, qq and qq* production for nig t q = 607 

GeV. The grey band shows the spread of unmatched Pythia predictions with varying shower 
parameters, while the green band shows the corresponding matched predictions. The full black 
curve represents the default unmatched Pythia Q 2 -ordered shower, while the full red curve is the 
default matched curve for gg production, included in the other plots as a guide for the eye. Each 
curve has normalization relative to the Born cross section, which in turn is normalized to unity in 
order to allow comparison between the different samples. 



show the total H? of radiated jets, denned as 

H^ = J2\P±i\ (3-1) 

where the sum is taken over radiated jets with p± > 40 GeV, as a measure of the total 
QCD radiation activity. 

An interesting observation from Fig. || is that there is no choice of shower parame- 
ters that describe all three types of particle production. For squark pair production in 
particular, the shower always undershoots the matched prediction, except for the hardest 
shower choice (j?j_-ordered power shower), which is much too hard in the other produc- 
tions. Indeed the matched curves are much more similar between the different productions 
than the shower curves, in particular in the tails (as illustrated by the inserted matched 
gg comparison curve). This indicates that the shower retains a too strong "memory" of 
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Figure 5: fl|? d and p 3 J_ distributions for radiated jets in gg production for m§ = 300,600, 1200 
GeV. The grey band shows the spread of unmatched Pythia predictions with varying shower 
parameters, while the green band shows the corresponding matched predictions. The full black 
curve represents the default unmatched Pythia Q 2 -ordered shower, while the full red curve is the 
default matched curve for rrig — 600 GeV, included in the other plots as a guide for the eye. Each 
curve has normalization relative to the Born cross section, which in turn is normalized to unity in 
order to allow comparison between the different samples. 

the initial state, while the matrix element displays a larger independence from the initial 
state. 

It is also interesting to note that the default virtuality-ordered Pythia shower under- 
shoots all curves, but least so for gluino production and most notably for squark-squark 
production, where it undershoots the matched curve by at least 50% starting already below 
a p± of 200 GeV. 

3.4 Produced particle mass dependence 

Determination of the absolute masses of produced particles at the LHC is suprisingly non- 
trivial if the produced particles decay to an invisible stable massive particle such as a 
WIMP. This difficulty is due to the fact that only the kinematics of decays are due to mass 
differences only, and are relatively independent of the absolute masses involved. Many 
methods have been devised to extract this information, including the use of different trans- 
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Scenarios 


I 


II 


III 


IV 




9 


600 


600 


heavy 


600 


Masses (GeV) 


Q 


550 


heavy 


550 


heavy 




X? 


100 


100 


100 


500 



Table 1: Benchmark scenarios employed in this work as modifications of SPSla 
assume the squarks decaying 100% into quark+lightest neutralino. 



We always 



verse mass quantities (see, e.g., [^, ^9|) and combined information on kinematics, cross 
sections and branching ratios [30, 31|, although especially the latter tend to be quite model 
dependent. Another possiblility would be to take advantage of the spectrum of QCD radia- 
tion, since this is directly sensitive to the absolute mass scale rather than mass differences. 
Studies have been done on the effect of QCD radiation in the decays of squarks and gluinos 
1 32, j33f| . In a hadron collider however, the majority of visible radiation originates from the 
initial state. 

As can be seen from Fig. ||, the predictions for the shapes of the spectra of radiated 
jets are quite precise once matching is taken into account, while it is clear that this type 
of study cannot be done with Pythia radiation only. 

For the lightest SUSY scenario (300 GeV gluinos), the distribution of the leading 
radiation jet is quite distinctive from the heavier scenarios, while the difference between 
the 600 GeV and 1200 GeV scenarios is small. Even barring the difficulties to distinguish 
radiation jets from jets from decays, it would be very difficult to differentiate between 
different high-mass scenarios using the distribution of radiated jets. One can however 
imagine scenarios where this could be a useful crosscheck, e.g. in the case where the mass 
splitting between the produced QCD particle and the stable WIMP is relatively small 



(such a scenario will be revisited in Sec. 44 below). More work on this would clearly be 
needed, but since this is beyond the scope for this paper, we choose to save it for a future 
publication. 



4. Anatomy of $± + multi-jet final states 

While the previous section dealt with jets from QCD radiation only, in models with new 
heavy QCD states produced at the LHC we expect jets also from their decays. It is therefore 
necessary to study the impact of jet matching in the context where we include also such jets. 
One could expect that jets from QCD radiation should be relatively unimportant compared 
to the hard decay jets from heavy particles, but as we will show, there are many situations 
when this is not true. Indeed, as it turns out, there are surprisingly many scenarios where 
radiation jets can significantly alter the analysis or interpretation of data. 

In order to clarify the discussion, and keep our results conservative, in this section we 
will use a set of simplified supersymmetric benchmark scenarios, summarized in Table |l|. 
In all the scenarios we assume all light-flavour squarks to have the same masses and that 
they all decay directly to the LSP, i.e. we ignore the existence of intermediate weak states. 
Introduction of cascade decays will have as main effect that jets from decays get softer, while 
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the jets from QCD radiation are not affected, and will hence mainly further accentuate our 
results. 

For ease of comparison between the scenarios, we have chosen to use the same masses, 
around 600 GeV, for the active heavy QCD states in all scenarios. Scenario I has a SUSY 
QCD spectrum similar to the SPS point la |p3|1 , with a gluino at 607 GeV which decays to 
squarks at 560 GeV, while the LSP is at 100 GeV. In scenario 2, the gluino has a mass of 
607 GeV but all squarks are heavy, so that the gluino decays through offshell squarks to 
two quarks and the LSP. Scenario III has squarks at 560 GeV and the gluino too heavy to 
be produced at the LHC, and finally scenario 4 has a gluino at 607 GeV decaying through 
offshell heavy squarks, but the LSP mass is 500 GeV, only 100 GeV lighter than the gluino. 
These scenarios will act as "cartoons" to illustrate different effects of QCD radiation in the 
production of new heavy QCD states. 

All plots and results in this section are generated using matched samples with the 
default Pythia parameter choices for virtuality-ordered showers unless otherwise stated. 

4.1 H 3 rp l variables in gluino production 

It is something of common lore that, in order to study squark pair production one should 
select 2-jet observables, for associated gluino-squark production at 3-jet observables and 
for gluino pair production at 4-jet observables. While it is obviously true that gluinos must 
decay to two quarks and a color singlet (barring the exotic possibility that the dominant 
decay is the loop-mediated twobody decay to a gluon and a color singlet), it should be kept 
in mind that the visibility of these jets depends strongly on the mass hierarchy of QCD 
states. In Fig. || we show the Hj! (n), defined as 

n 

4*(n) = X>if I (4.1) 

i=l 



for n = 2,3,4, for the scenarios 1 and 2. The sum in eq. (4.1) is taken over jets defined 
using the SISCone |l5| algorithm with a radius of 0.5 and p± > 40 GeV. We also show the 
composition of (n) in terms of jets from the gluino decay and radiated jets ("ISR"), as 
well as show the average fraction of the coming from the decay. For scenario 2, where 
the gluinos decay through off-shell squarks to two quarks of similar energies and an LSP, 
the majority of events in the peak of i7^ t (4) include only jets from the decay, while the 
tail of the distribution is dominated by 2-3 jets from the decay and one jet from radiation. 
For scenario 1 however, where the 600 GeV gluino decays into a fairly soft jet (with an 
energy around 50 GeV) and a squark which in turn decays to a hard jet and an LSP, the 
distributions are quite different. Here, ff^? t (4) is dominated by events where at least one 
of the jets come from radiation, and we need to go down to H^ l (2) to be dominated by 
events with only decay jets across the whole Ht range. The energy fraction of iTjT (4) 
coming from decay is still fairly high, even in the high-ify tail, since most of the transverse 
energy comes from the squark decay jets. 

The immediate interpretation of this result is that in a scenario with a small mass 
splitting between gluinos and squarks, gluino production might be difficult to distinguish 
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Figure 6: (n) for n — 3,4, 5 in gg production. The different point markers show how many 
of the jets entering in the definition of iJ^'(n) come from QCD radiation. Upper row: Scenario I, 
lower row: Scenario II. Below each main plot, the red curve indicates the percentage of Ht coming 
from the decays. Jets are defined using the SISCone algorithm with a p™ m of 40 GeV and a radius 
of 0.5. 

from squark production with additional QCD radiation. In Section |4.3[ we will study this 
in greater detail. 

4.2 Jet multiplicities for different scenarios 

A question that the proper matching of jets is particularly apt to answer, is to specify the 
number of jets typically present in production of different particles. We here present a table 
with the jet multiplicities, for matched and unmatched (Pythia virtuality-ordered default) 
production, for the scenarios studied. In order to make the table as useful as possible, we 
have used jet cuts close to what is used in many preparatory analyses for squark and gluino 
searches: pf 1 > 180 GeV, pf 2 > 110 GeV, p 3 ± j > 50 GeV for j > 2 and \r]i ct \ < 3 for all 
jets. We require all events to have at least two jets. The jet multiplicities are exclusive, and 
so add up to 100% of the events passing the 2-jet cut. The jet algorithm used is SISCone 
with a radius of 0.5 and p± > 40 GeV. 

The difference between matched and unmatched generation (as seen in Fig. is 
twofold: First matching tends to increase the transverse boost of the produced pair, and 
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Process 


AT 
1 V 


— Z 


AT 
iv 


_ Q 
— O 


AT 


— A 

— 4 


AT 
1 V 


— K 

— 


N ^ 6 


Signal eff. 




M 

1V1 


TT 

U 


M 


TT 

U 


1V1 


TT 

U 


M 


TT 

U 


M 


U 


M 


U 


99 sc. I 


15.7 


27.1 


30.0 


33.9 


24.4 


21.0 


13.0 


8.9 


16.7 


8.9 


43.7 


40.4 


gq sc.I 


35.2 


39.4 


32.5 


33.8 


17.5 


16.0 


7.0 


5.2 


7.64 


5.2 


31.9 


28.3 


qq sc.I 


40.2 


48.1 


33.4 


32.0 


15.9 


12.9 


4.9 


3.8 


5.4 


3.4 


16.9 


16.0 


99 sell 


4.4 


4.7 


19.5 


22.1 


27.1 


29.2 


18.0 


17.6 


31.1 


26.1 


43.8 


40.1 


99 sc.IV 


21.5 


28.4 


32.6 


37.0 


23.9 


21.0 


11.4 


7.8 


10.4 


5.8 


4.7 


3.0 



Table 2: Contribution of events with N jets for matched and unmatched processes: gg, gq and qq 
in Scenario I (for squark production this is very similar to scenario III), and gg in Scenario II and 
IV. All numbers are in percent. "Signal efficiency" shows the percentage of events that pass the 
2-jet cut. The jet cuts are described in the text. 



hence the p± of the softest jets from the decays. Second, the ISR jets get harder, and more 
easily get above the threshold of 50 GeV. The matched production therefore in general 
populates higher jet number bins the unmatched ones. 

The first three lines in Table ^ represent Scenario I, which is very similar to the 
benchmark point SPSla. In this scenario, the gluinos decay as g — > qq — ► qqx° (with 
m x o ~100 GeV), so there are typically two hard and two soft jets from the gluino decay. 
Here the effect of the matching is large, due to the increase in p± for the soft jets as well as 
increased hardness of the ISR jets. For scenario II on the other hand, where gluinos decay 
in three body decay into two jets and a light neutralino, the sensitivity to the matching is 
much lower, since there, typically at least four reasonably hard jets are present from the 
decays. It is only in this and similar cases that the statement that gluino pair production 
generally corresponds to four hard jets in the event is true. In the last row in the table, 
the produced gluinos decay to a near-degenerate LSP, meaning that most hard jets are due 
to QCD radiation. Here, a large recoil against initial state jets is needed in order to even 
pass the 2-jet cut, hence the very low signal efficiency. 

For squark-squark and gluino-squark production for Scenario I (row 2 and 3 in Table [2], 
there are only two hard partons from the decays (and one additional soft parton from the 
gluino decay) and additional QCD radiation. Again the main consequence of the matching 
is to increase the mean number of jets. It is interesting to note that in this particular case, 
the addition of matching to the generation of squark pair production gives very similar 
numbers as the unmatched gq associated production (although the selection efficiency is 
different), indicating that the matching has a similar impact as the addition of one extra 
jet with p± ~ 50 GeV. 

We now look closer at two examples where an analysis based only on a parton shower 
approach might lead to the wrong conclusions. 

4.3 Example 1: False gluino evidence 

We here consider the case given in Scenario III, where the only observable SUSY particles 
are 560 GeV squarks, while the gluinos are too heavy to be produced. If such a scenario is 
simulated using only parton showers, the result will typically be a deficit in the number of 
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Figure 7: Ht{ti) with n — 2,3,4 for qq^*' production in Scenario III, with and without matching 
(full black and full red curves), and for unmatched associated gq and gg production with nig — 700 
GeV(dotted and broken red curves). Also shown is the result of adding the unmatched qq(*' curve 
and 0.25 x the gq curve, which perfectly mimicks the matched curve. All simulations are done with 
Pythia default Q 2 -ordered showers. 

events with multiple jets. This is illustrated in Fig. [7], where we show (2), H^ t (3) and 
H^(A) (defined in eq. ( |4.lD ) for matched and unmatched simulation of qq^ production, 
using the Pythia virtuality-ordered shower with default parameters (the full black and 
full red curves, respectively). While the unmatched generation well reproduces the 2-jet 
Ht, it increasingly falls below the matched curve for H^{3) and 11^(4). In real data, 
this deficit could easily be interpreted as a sign for a missing production mode, such as 
associated production with gluinos and squarks, or gluino pair production. 

Also shown in Fig. [7| are the curves for associated gluino production (dotted) and gluino 
pair production (broken red), with gluinos at 700 GeV, decaying to the squarks emitting an 
additional jet. While the cross section for the associated gq production process is too large 
to be accommodated by the "data" (i.e. the matched squark production curve), there are 
many ways in which this cross section estimate could be wrong; if, for example, the decay 
modes are different for squarks and gluinos, or if the particles seen are not from the MSSM 
but from some other realization of new physics. We can therefore view the normalization 
of the different curves as free parameters, and look at how the curves could be combined 
to generate a fit to the "data" curve, i.e. the matched qcfi*' curve. 

It turnes out that for these choices of masses, the data can be very well accommodated 
by adding the gluino-squark associated production, with the cross section reduced to 25% 
of the nominal cross section, to the unmatched squark production curve, as shown by the 
broken blue line in the figure. 

Possible cross checks to avoid this type of false conclusions could be looking at detailed 
kinematic properties of the jets, or looking at different decay signatures. However, this 
example still illustrates the importance of using matching also for BSM signals, in particular 
when looking at observables where extra jets from QCD radiation may play an important 
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Figure 8: ft±_ and p± for the two hardest jets in gg production in Scenario IV. The grey band 
shows the spread of unmatched Pythia predictions with varying shower parameters, while the 
green band shows the corresponding matched predictions. The full blue and red curves represents 
the default unmatched Pythia Q 2 - and pj_ -ordered showers respectively. 



role. 

4.4 Example 2: Degenerate spectrum 

Another case where jet matching is of considerable importance is when the mass splitting 
between the produced QCD particle and the LSP is small, as is the case in Scenario IV. 
There, 600 GeV gluinos decay through off-shell squarks to a neutral 550 GeV LSP and 
two quarks. The main problem with searches in this scenario is that no large missing Et 
is produced in the decay, since all decay products are soft. It doesn't help if the gluinos 
are produced at a large invariant mass, since they will then decay back-to-back, and the 
missing energy due to the boost of the individual gluinos is canceled between the two LSP's. 
While this type of scenario is not possible in typical unification scenarios (e.g. mSUGRA 
or mGMSB), a more model-independent approach necessitates taking them into account 

MM- 

The only way to get a large missing energy in this scenario is through a transverse boost 
of the gluino pair center-of-mass system. This happens when the gluino pair system recoils 
against hard initial state radiation. In this case, a proper jet matching between parton 
showers and matrix element QCD emission is crucial to well describe the E± distribution 
as well as the resulting jet structure. This is particularly important since a E± of a certain 
magnitude as well as several hard jets is typically needed to pass the experiment triggers, so 
that a misrepresentation of these quantities will lead to very different trigger acceptences. 

We show in Fig. |8| three representative quantities for this case, the B± and the p± of 
the two hardest jets. As can be seen from the figure, the virtuality-ordered Pythia default 
distribution severely underestimates the missing Et as well as the hardness of the leading 
jets, which mainly come from the initial state radiation. The effect, with typical cuts jet- 
E± cuts like E± > 200 GeV, pf 1 > 180 GeV and pf 2 > 110 GeV is a signal efficiency less 
then half that of the matched production. The p^-ordered shower rather overshoots the 
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matrix element curves, up to the factorization scale (m g ), where it falls off rapidly. Once 
again, the lesson here is that in order to get a description that is predictive and insensitive 
to the details of the shower parameterization, it is necessary to use jet matching. 



5. Impact on BSM searches 

We now consider the effects of an accurate simulation of QCD radiation in the typical 
observables employed in the BSM inclusive searches: high p± jets and high missing trans- 
verse energy. In order to see the effect of the matching on the sensitivity to shower after a 
smearing of the signal, and in order to be as complete as possible, we consider the produc- 
tion of gg, gq, qq(*> and tx^t* 2> with q defined as ul,r, di,R, sl,r, cl,r, 61,2- The signal is 
produced in both the matched (2— >2, 3, 4) and unmatched (2— >2) modes in the the SPSla 
benchmark scenario [34]. For the background, we consider the most important processes 



leading to four hard jets and potentially large missing transverse energy: — > l^u + 4 
jets, Z° — > vv + 4 jets, — > Tj e tv + 3 jets and finally the inclusive ti+0,1,2,3 jets. We do 
not include QCD multijet production, since we have no means of realistically performing 
simulations of the missing energy distribution, which is due to decays of heavy quarks 
to neutrinos and jet mismeasurement in the detector. We instead base our analysis on 



cuts similar to those used in Refs. 38] and [39], and keep this contribution in mind. All 
background simulations are done using jet matching. Many comparisons have been done 
between matched and unmatched background simulations, which are well known to differ 
by up to several orders of magnitude for this type of multi-jet observables JTJ, |(], |T^| . We 
therefore here look at the effects of including matching only in the signal simulation. 

Detector simulation is performed using PGS 4 H with the MidPoint cone algorithm 
with a minimum p± of 40 GeV and a radius of 0.5. To be conservative, we use only 
kinematic variables associated to the jets and the missing transverse energy. The cuts used 
are 

• A/jet > 4 

• |t7i| < 1.7, [7/2,3,4... I < 3 

• pf 1 >180 GeV, p^ 2 > 110 GeV, pf >2 > 50 GeV 

• #L > 150 GeV 

• Acf)($ ± ,}et 1 ) > 0.5 and A^(^j_,jet 2 ) > 1 

• Et=2Pf l +^± > 600 GeV . 

The SPSla scenario is affected by several of the difficulties described in Sec. ^. The 
gluino has a mass higher than, but close to, the squarks, and hence decays to a squark 
and a soft jet, making the jet counting complicated. Since gluinos and squarks are of 
similar mass, the QCD SUSY production includes associated squark-gluino production, 
gluino pair production and squark pair production (both qq and qg* which are of similar 
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Figure 9: Ht for SUSY production in the SPSla benchmark scenario, compared to (matched) 
Standard Model backgrounds. The grey band shows the spread of unmatched Pythia predictions 
with varying shower parameters, while the green band shows the corresponding matched predictions. 
Ht definition and cuts are described in the text. 

cross section, and to a negligible degree q*q*), in order of cross section. This means that 
the Pythia shower cannot simultaneously describe all production modes, as demonstrated 
in Fig. f|. Furthermore the separation of the different production modes based on jet 
counts or jet kinematics will be non-trivial. The production cross section is however in 
this scenario dominated by associated gluino-squark and gluino pair production, where 
the default Pythia description is reasonably good, and only undershoots the matched 
description by about 10-30% for the first couple of QCD radiation jets, so we expect the 
inclusive Pythia description to be reasonably close to the matched curve. 

The result is illustrated in Fig. ^, which shows the Ht = X^=2 Pr ' + $± Ior inclusive 
supersymmetric production of gluinos and squarks (including t and b). 

The effect of the matching is as expected - a significant reduction in the sensitiv- 
ity to parton shower parameters, and a shift of the prediction as compared to the default 
virtuality-ordered Pythia shower by about 10-30%. Even with the smearing due to the de- 
tector simulation and more complex decays than the simplified scenarios used in Sec. the 
"power" pj_-ordered shower continues to overshoot the result obtained with the matching 
whereas the "wimpy" virtuality-ordered showers undershoot the matched curve. 

Since we require four hard jets, the strongest impact is on the squark-(anti)squark 
pair production, with an efficiency increase close to 40% when passing from the unmatched 
default virtuality-ordered Pythia shower to the matched production. This happens for 
two reasons; as described in Sec. [h^, only two hard jets are produced by the decay of the 
squarks, which means that two jets from QCD radiation are needed. Second, as shown in 
Fig. |], the difference between the unmatched and matched radiation is particularly large for 
squark production. On the other hand, the effect is much smaller for gluino pair production 
(around 5%), which has to do mainly with the large fraction of events with at least one 
top quark in the decay, giving rise to multiple hard jets and hence a small sensitivity to 
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the matching. The effect on the associated gluino-squark production lies between these 
two extremes, with an efficiency increase close to 10%. t and b production is dominated 
by tit* pairs (due to a low i\ mass, around 390 GeV), which decay to top quarks and the 
LSP. For those, the effect of the matching is negligible. 



6. Conclusions 



Discovering new physics at the LHC will most probably be quite challenging. Apart from a 
few very clean but theoretically not very motivated signatures, such as narrow resonances 
decaying into lepton pairs, the vast majority of the theoretical constructions which are still 
viable nowadays gives a rich and complex phenomenology that will be hard to decipher. For 
instance, many models that address the hierarchy problem, predict new heavy colored states 
at the TeV scale that can be easily produced at hadron colliders. To what extent the decay 
products of such states will be identifiable on top of a large QCD background, including 
vector boson(s) + jets, ti jets and multi-jets, is currently subject of many theoretical and 
experimental investigations. 

In this work we have for the first time addressed the issue of quantifying the effects of 
extra QCD radiation in the production of heavy colored states employing inclusive multi-jet 
samples obtained by matching multiparton matrix elements and parton showers. Previous 
work has either considered radiation in Standard Model tt events only or in SUSY but 



at parton level [25|. 

Our results can be briefly summarized as follows. First the extension of matching 
techniques to beyond the Standard Model scenarios, such as SUSY or UED, while posing 
no problems of principle, requires dealing with new technical issues. For instance, the 
combination of multi-parton samples for production of different resonances decaying into 
jets leads to problems of double counting that need to be addressed. We have proposed 
two working solutions, both of which have been implemented in MadGraph/MadEvent. 

The main result of this work is that a matched matrix element plus parton shower 
approach for heavy particle production is in general much more accurate and predictive 
than a parton shower alone. We find that, contrary to the "common lore" that showers 
alone provide a good description of extra radiation when the produced states are heavy, 
there are many cases where matrix element corrections are indispensable. Not only are 
the p± spectra of the extra jets in the parton shower approach extremely sensitive to the 
shower starting scale and the shower algorithm, while the matched simulations are not, 
but also that there is in general no tuning of the shower that can simultaneously reproduce 
the matched samples for all initial states. We have presented examples where simulations 
based only on the parton shower could be misleading. Another possibly relevant issue, 
that we leave to future studies, is matching of radiation in the decays of the heavy states. 
Though already possible in our current implementation, we have not included it in this 
study mainly because it is computationally expensive. Next-to-leading order effects in 
decays have previously been shown to be possibly important in precision studies, e.g. of 
the spins of new particles |33|]. If found relevant, also matching of radiation in decays could 
be included in future simulations. 
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In conclusion, we recommend that matched samples should be used not only for back- 
grounds but also for beyond the Standard Model physics signals. Having the best available 
simulations could be important not only in designing better and more solid strategies to 
make discoveries but even more to identify the nature of the new physics. For precision 
measurements of new physics properties, as well as reliably distinguish between different 
scenarios, fully differential next-to-leading order simulations will most probably be needed. 
Given the large number of theoretical possibilities still open, it is clear that an automatic 
approach to NLO computations, and to their matching with a parton shower, will be 
certainly welcome. 
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